function vb = glm_VOIBetas(hfile, vfile, opts)
% GLM::VOIBetas  - returns a table of VOI betas (per subjects)
%
% FORMAT:       vb = glm.VOIBetas(voi [, opts]);
%
% Input fields:
%
%       voi         VOI object
%       opts        optional settings
%        .c         contrasts (defaults: each beta map on its own)
%        .interp    either of {'nearest'}, 'linear', 'cubic'
%        .rmean     remove mean (of map) first
%        .robust    perform robust mean estimation of average (for VOIs
%                   with at least 4 voxels)
%        .vl        indices of VOIs in object to sample (default: all)
%
% Output fields:
%
%       vb          SxCxV double table with data

% Version:  v0.7f
% Build:    8110721
% Date:     Nov-07 2008, 9:11 PM CET
% Author:   Jochen Weber, SCAN Unit, Columbia University, NYC, NY, USA
% URL/Info: http://wiki.brainvoyager.com/BVQXtools

% check arguments
if nargin < 2 || ...
    numel(hfile) ~= 1 || ...
   ~isBVQXfile(hfile, 'glm') || ...
    numel(vfile) ~= 1 || ...
   ~isBVQXfile(vfile, 'voi')
    error( ...
        'BVQXfile:BadArgument', ...
        'Invalid object handle in call.' ...
    );
end
bc = bvqxfile_getcont(hfile.L);
if bc.ProjectType ~= 1
    error( ...
        'BVQXfile:BadArgument', ...
        'Only valid for VTC based GLM files.' ...
    );
end
switch (bc.ProjectTypeRFX)
    case {0}
        ns = 1;
        np = bc.NrOfPredictors;
    case {1}
        ns = numel(bc.GLMData.Subject);
        np = bc.NrOfSubjectPredictors;
    otherwise
        error( ...
            'BVQXfile:BadArgument', ...
            'Invalid/unsupported ProjectTypeRFX flag in file.' ...
        );
end
vc = bvqxfile_getcont(vfile.L);
if nargin < 3 || ...
   ~isstruct(opts) || ...
    numel(opts) ~= 1
    opts = struct;
end
if ~isfield(opts, 'c') || ...
   ~isa(opts.c, 'double') || ...
    isempty(opts.c) || ...
   (size(opts.c, 2) ~= np && ...
    size(opts.c, 2) ~= (np - 1)) || ...
    any(isinf(opts.c(:)) | isnan(opts.c(:))) || ...
    any(any(opts.c < 0, 2) & sum(opts.c, 2) ~= 0)
    opts.c = eye(np);
end
nc = size(opts.c, 1);
if ~isfield(opts, 'interp') || ...
   ~ischar(opts.interp) || ...
    isempty(opts.interp) || ...
   ~any(lower(opts.interp(1)) == 'cln')
    opts.interp = 'n';
else
    opts.interp = lower(opts.interp(1));
end
if ~isfield(opts, 'rmean') || ...
   ~islogical(opts.rmean) || ...
    numel(opts.rmean) ~= 1
    opts.rmean = false;
end
if ~isfield(opts, 'robust') || ...
   ~islogical(opts.robust) || ...
    numel(opts.robust) ~= 1
    opts.robust = false;
end
if ~isfield(opts, 'vl') || ...
   ~isa(opts.vl, 'double') || ...
    isempty(opts.vl) || ...
    any(isinf(opts.vl(:)) | isnan(opts.vl(:)))
    opts.vl = 1:numel(vc.VOI);
else
    opts.vl = intersect(opts.vl(:)', 1:numel(vc.VOI));
end

% prepare output
nv = numel(opts.vl);
vb = zeros(ns, nc, nv);

% get bounding box
bb = aft_BoundingBox(hfile);

% convert coordinates
voi = cell(1, numel(opts.vl));
if strcmpi(vc.ReferenceSpace, 'tal')
    convtype = 'tal2bvc';
else
    convtype = 'bvs2bvc';
end
for c = 1:numel(voi)
    voi{c} = bvcoordconv(vc.VOI(opts.vl(c)).Voxels, convtype, bb);
    if opts.interp == 'n'
        voi{c} = round(voi{c});
    end
end

% interpolation option
switch (opts.interp)
    case {'c'}
        ipo = 'cubic';
    case {'l'}
        ipo = 'linear';
    case {'n'}
        ipo = 'nearest';
end

% for RFX GLMs
if bc.ProjectTypeRFX == 1
    
    % iterate over subjects
    for sc = 1:ns
        
        % iterate over contrasts, predictors and VOIs
        for cc = 1:nc
            for pc = 1:size(opts.c, 2)
                cval = opts.c(cc, pc);
                if cval == 0
                    continue;
                end
                bm = bc.GLMData.Subject(sc).BetaMaps(:, :, :, pc);
                bn = (bm ~= 0);
                bm(~bn) = Inf;
                if opts.rmean && ...
                    pc < np
                    bm(bn) = bm(bn) - mean(bm(bn));
                end
                for vc = 1:nv
                    if ipo(1) ~= 'n'
                        ipv = flexinterpn_method(bm, voi{vc}, ipo);
                    else
                        ipv = bm(sub2ind(size(bm), ...
                            voi{vc}(:, 1), voi{vc}(:, 2), voi{vc}(:,3)));
                    end
                    ipv(isinf(ipv) | isnan(ipv) | ipv == 0) = [];
                    if ~opts.robust || ...
                        numel(ipv) < 4
                        vb(sc, cc, vc) = vb(sc, cc, vc) + cval * mean(ipv);
                    else
                        vb(sc, cc, vc) = vb(sc, cc, vc) + cval * ...
                            fitrobustbisquare(ones(numel(ipv), 1), double(ipv));
                    end
                end
            end
        end
    end
else
end
